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Abstract.  We  have  focused  upon  the  development  and  validation  of  finite 
element  methods  for  LES  of  turbulent  flows  in  settings  in  which  interac- 
tion with  (possibly  geometrically  complex)  boundaries  are  important.  New 
results  are  presented  in  Section  2 on  closure  and  convolution  on  bounded 
domains.  Similarly,  new  near  wall  models  were  required;  our  method  for 
developing  these  is  described  and  one  from  our  work  is  presented  in  Section 
3.  The  difficulties  in  closure  and  wall  modelling  suggest  a second  approach: 
direct  simulation  of  large  eddy  motion.  We  give  an  extension  of  this  ap- 
proach to  nonlinear,  equilibrium  flows  - a step  closer  to  turbulence  for  an 
approach  not  requiring  wall  models  or  closure  models. 


1.  Introduction 

There  is  a natural  interplay  between  Finite  Element  CFD  and  LES  for  top- 
ics such  as:  closure  models,  near  wall  models,  and  FEM  postprocessing  by 
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local  averaging.  Because  of  their  geometric  flexibility  and  their  flexibility 
with  respect  to  continuum  models  and  linear  or  nonlinear  boundary  condi- 
tions, FEM’s  are  a natural  discretization  in  the  LES  of  flows  with  complex 
boundaries.  Adaptivity,  highly  developed  for  FEM’s,  has  interesting  po- 
tential for  LES  without  modelling  error.  For  example,  using  extensions  of 
work  in  (John  and  Layton,  2001)  reported  in  Section  4,  a CFD  mesh  can 
be  designed  so  that  local  averages  of  an  inaccurate  approximate  flow  field 
approximate  the  true  local  velocity  averages  with  assured  accuracy. 

This  paper  surveys  some  of  our  work  developing  a mathematical  foun- 
dation for  finite  element  LES  and  presents  some  new  results  and  extensions 
of  this  work.  Many  of  the  reports  cited  herein  are  available  at 
http://www.math.pitt.edu/  rsj  wjl. 

Consider  therefore  the  turbulent  flow  of  an  incompressible  fluid  in  three 
dimensions,  bounded  by  walls  and  driven  by  a body  force  f{x,t).  The 
velocity-pressure  (u,p)  satisfy  the  Navier  Stokes  equations,  given  by: 

ut  + div(uu)  — 2iie_1div  (B(u))  + gradp  = f(x,  t),  and  divu  = 0,  (1) 

in  a domain  Q complemented  by  boundary  and  initial  conditions 
n = 0on  walls,  and  u(x,  0)  = uq{x), 

where  O(u)  = (grad  u+  grad  ut)/2  is  the  velocity  deformation  tensor.  Pick- 
ing a length  scale  6 and  an  associated  averaging  kernel  gs{x)  :=  6~3g(x/6), 
where  g(x)  is  a mollifier  satisfying  certain  properties,  local  velocity  and 
pressure  averages  are  frequently  defined  by  convolution  with  g{x),u  := 
gs  * u,p  :=  gg  *p,  etc.,  where  all  functions  are,  when  necessary,  extended 
by  0 off  the  flow  domain  to  compute  the  required  average. 

For  compactness,  we  will  only  discuss  herein  the  case  of  constant  aver- 
aging radius  S.  Alternate  approaches  include  a variational  definition  of  the 
large  eddies  (Layton,  1999),  (Hughes  et  al.,  2001),  and  differential  filters, 
introduced  by  Germano  in  the  1980’s,  wherein  (in  effect)  <5  = 6(x)  — » 0 
as  x approaches  walls  in  a manner  intrinsic  to  the  NSE  and  the  domain’s 
geometry,  (Layton  and  Lewandowski,  2001). 

2.  Convolution  on  Bounded  Domains  and  Closure 

With  constant  averaging  radius,  6,  convolution  operators  commute  with 
differential  operators  in  the  absence  of  boundaries.  With  boundaries,  extra 
terms  arise  which  are  often  overlooked.  For  example,  filtering  (1)  with  gs(x) 
on  a bounded  domain  gives  the  following  space  filtered  equations  for  u := 
gs  * u,p  :=  gs  * p (see  (Dunca  et  al.,  2001)) 

ut  + div(tt  u)  — 2/?e-1div(B(tt))  + gradp  + div(T(u))  = / + As{o),  (2) 
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where  o is  the  stress  of  the  unknown  flow  and  T(u)  :=  uu  — u u is  the 
Reynolds  stress  tensor.  A direct  calculation  in  (Dunca  et  al.,  2001)  using 
the  theory  of  distributions  reveals  that 

As  (a)  :=  / gs(x  — s)(a  ■ h)(s)ds,  a :=  2Re~lH(u)  — pi, 
an 

which  can  be  interpreted  as  a boundary  stress  distribution  smeared  out  onto 
the  entire  flow  domain.  Since  the  As  (a)  term  in  (2)  is  normally  neglected, 
it  might  be  hoped  to  be  negligible.  Unfortunately,  this  is  not  the  case,  as  is 
seen  in  the  following  result,  (Dunca  et  al.,  2001). 

Proposition.  Let  ||m||x,P  :=  (f  \w\pdx)l!p , 1 < p < oo  denote  the  usual  Lp 

R3 

norm,  1 < p < oo  and  ||itf||L°°  :=  ess  sup  |m|.  Then  for  1 < p < oo 

zeR3 

H-^llz/P  ~ t 0 as  6 — y 0. 

if  and  only  if  the  normal  stress  is  identically  zero  on  the  boundary  of  the 
flow  domain 

a ■ n = 0 on  the  boundary.  □ 


The  proof  of  this  result  is  moderately  technical  but  the  importance 
of  the  result  is  clear:  the  term  A$(cr)  cannot  be  omitted  if  the  boundary 
influences  the  flow! 

It  is  well-known  that  (2)  are  not  closed  for  periodic  problems  due  (only) 
to  the  Reynolds  stress  tensor  T (u)  and  numerous  LES  models  have  been  de- 
veloped to  model  T(tt),  see,  e.g.,  (Iliescu  and  Layton,  1998),  (Berselli  et  al., 
2001),  (Layton,  2000),  (Galdi  and  Layton,  1999),  (Layton  and  Lewandowski, 
2001)  for  our  work  on  closure.  For  flows  with  real  walls  or  other  boundaries, 
(2)  are  not  closed  due  to  both  T(u)  and  the  “smeared”  wall  stress  oper- 
ator Ag(a).  More  detailed  understanding  of  As(a)  is  necessary  to  model 
correctly  the  interaction  of  large  eddies  with  walls;  see  (Dunca  et  al.,  2001) 
for  first  steps. 

3.  Near  Wall  Models:  Boundary  Conditions  for  LES 

Either  modeling  or  omitting  the  smeared  wall  stress  term  Aj(-),  boundary 
conditions  are  still  required  for  the  large  eddies  u and  simply  imposing 
u = 0 at  walls  is  inconsistent,  see  (Galdi  and  Layton,  1999),  Figure  1.  This 
inconsistency  is  intuitively  clear:  a tornado,  as  an  example  of  a large  eddy, 
does  inove/slip  along  the  ground. 
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What  are  then  the  correct  boundary  conditions?  Motivated  by  this  ex- 
ample of  a tornado,  we  have  explored  no-penetration  and  slip-with-friction 
conditions  (see,  e.g.,  (Navier,  1823)  and  (Maxwell,  1879)  for  antecedents) 
for  u: 

u • h = 0 and  (3 j u-  fj  + t(u)  ■ fj  = 0 on  walls,  (3) 

where  n, ti,T2  are  unit  normal  and  tangent  vectors  to  the  wall  (Layton, 
2001),  and  t is  the  Cauchy  stress  or  traction  vector.  The  effective  friction 
coefficient  (3  for  the  large  eddies  is  calculated  explicitly  in  (Sahin,  2000) 
and  (John  et  al.,  2001) 

(3  = (3(8,  Re)  or  (3  = (3(8,  | u ■ fj\)  : large  eddy  friction  coefficient. 

The  conditions  (3)  thus  allow  both  linear  and  nonlinear  near  wall  models. 
For  turbulent  channel  flow,  the  linear  model,  based  on  a global  Reynolds 
number,  seems  to  suffice.  However,  for  flows  in  complex  geometries,  the 
local  Reynolds  number,  related  to  the  slip  velocity  \u  ■ fj\,  varies  greatly 
from  recirculation  zones  to  mean  stream  regions.  Thus,  the  nonlinear  model 
seems  necessary. 

The  friction  laws  in  (Sahin,  2000),  (John  et  al.,  2001)  are  calculated 
using  boundary  layer  theory,  with  different  results  for  different  types  of 
turbulent  layers.  For  example,  for  a power  law  layer  the  linear  wall  model’s 
effective  friction  coefficient  is  as  follows.  Let  T[-]  denote  the  usual  gamma 
function  and  T[-,  •]  the  incomplete  gamma  function,  rj  = 0.21  Re-1  and 
Z:=r)/8 


a,s  „>  = Jfe--z'/7[7r[jj]-ir[ji,z3]l 

l,W  2vW[Z1'7{r[f]-r[f,Z2]}  + V?(l  - erf{Z))] 

A very  good  (and  simple)  approximation  is  given  (to  4 decimal  places)  by 
the  effective  friction  coefficient: 

(3(8,  Re)  “ i?e-1<5-21.22e~°'008/z.  (4) 


From  this  simple  formula  we  see  that: 

— As  8 — > 0,  for  fixed  Re,  the  friction  (3  — > oo,  i.e.,  the  boundary  condi- 
tions (3)  reduce  to  no  slip. 

— As  Re  -»  oo  for  fixed  8,  the  friction  coefficient  (3  — > 0,  i.e.,  the  boundary 
conditions  (3)  reduce  to  free-slip  and  no-penetration,  as  appropriate  for 
the  Euler  equations. 

To  explain  the  nonlinear  friction  laws  we  note  that  for  a power  law 
layer  (and  other  profiles  as  well),  the  slip  speed,  s :=  y/\ u • f\\2  + \u  ■ faj2 
is  a monotone  increasing  function  of  Re  for  8 fixed.  Expressing  this  as 
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5 = g{Re),g'(Re)  > 0 implies  that  the  inverse  function  <?*(•)  of  g(-)  exists. 
We  can  then  write  Re  = g*(s).  Inserting  this  inverted  relation  into  (4) 
yields  nonlinear  friction  laws: 

Pj{6,g*{\u  ■ Tj\))  “ 1.22 rV(|u  • fjDe"0’0090/®*!1®^1). 

This  program  can  be  carried  out  for  different  types  of  turbulent  boundary 
layer  profiles,  see  (John  et  al.,  2001)  for  detailed  formulas. 

4.  Numerical  Errors  in  LES 

With  closure  models  and  boundary  conditions  selected,  convergence  of  a 
chosen  numerical  method  to  the  model’s  solution  must  be  considered,  e.g. 
(John  and  Layton,  2000),  (Iliescu  et  al.,  2000),  (Layton,  1996).  For  example, 
in  (John  and  Layton,  2000)  convergence  for  a FEM  discretization  of  the 
Smagorinsky  model  is  proven  to  be  uniform  in  Re,  as  is  often  reported  in 
practical  simulations  but  hitherto  unproven.  This  work  does  not  address 
the  modelling  error  itself  however. 

In  recent  work,  we  have  devised  new  algorithms  which  give  computable 
bounds  on  both  the  modelling  and  the  numerical  errors,  (John  and  Layton, 
2000).  This  work  was  begun  in  (John  and  Layton,  2001)  for  (linear)  Stokes 
flow.  Our  current  work,  reported  in  this  section,  has  extended  those  new 
methods  to  the  (nonlinear)  equilibrium  NSE  - a step  closer  to  the  true 
problem. 

The  approach  is  simple  to  describe:  the  fluid  velocity  is  directly  ap- 
proximated and  then  postprocessed  by  local  averaging,  giving  g$  * uh  as  an 
approximation  to  u — g$  * u.  One  key  is  that  the  finite  element  mesh  is 
adapted  within  the  calculation  so  that  gg  * uh  approximates  u with  assured 
accuracy,  even  when  uh  is  a bad  approximation  to  u.  The  second  key  is  that 
the  averaging  radius  S must  be  taken  smaller  than  the  local  meshwidth  of 
the  refined  mesh,  6 « O(h).  These  new  a posteriori  estimators  are  given 
in  (John  and  Layton,  2000)  (and  a new  report  in  preparation).  With  these 
new  estimators,  the  error  is  also  concentrated  in  high  frequencies  and  killed 
by  postprocessing.  Thus,  the  error  in  the  large  eddies  is  typically  far  smaller 
than  the  overall  error,  as  in  the  following  theorem. 

Theorem.  Let  the  equilibrium  NSE  be  solved  by  the  FEM  with  finite  ele- 
ment spaces  satisfying  the  usual  stability  and  local-polynomial  degree  k ap- 
proximation conditions.  Let  the  NSE  solution  be  nonsingular  with  linearized 
dual  Hk+1 -regular.  Then,  the  error  in  the  large  eddies  \\gs  * uh  — u||p  is 
related  to  ( and  much  smaller  than),  the  error  in  the  velocity  (||V(it  — u"')||) 
and  pressure  (||p  — ph\\)  by,  for  any  e > 0: 

\\gs  *uh -u\\L2  < + y/hj  (^j  <f||V(u  - uh)\\L2 
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+ (jj  ■5(l|p-P*l|t>  + l|V-(«-«',)IU.) 

+ (f  «I|V(»  - uh)\\l,  + «H||„  - □ 

This  theorem  also  gives  analytical  guidance  for  relating  S and  h.  How- 
ever, the  related  mesh-adaptation  strategy  is  the  most  important  practical 
contribution  of  this  theory. 
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